##############################
#Define functions for looping#
##############################

#Loop descriptive statistics
summ_loop = function(y, data) {
  x = data[[y]]
  data.frame(
    Variable = y,
    min = min(x, na.rm = T),
    max = max(x, na.rm = T),
    mean = mean(x, na.rm = T),
    sd = sd(x, na.rm = T),
    N = x %>% na.omit %>% length
  ) %>% return
}

#Loop least squares regression
lm_loop = function(y, x, data, w = NULL){
  f = as.formula(paste(y,"~",x, collapse = ""))
  out = lm(f, data = data, weights = w)
  return(out)
}

#Loop negative binomial regression
glmnb_loop = function(y, x, data){
  require(MASS)
  f = as.formula(paste(y,"~",x, collapse = ""))
  out = glm.nb(f, data = data)
  return(out)
}

#Clustered standard errors
cluster_errors = function(model, clusters, vcov = T) {
  require(multiwayvcov)
  require(lmtest)
  cl_vcov = cluster.vcov(model, clusters)
  
  if (!vcov) {
    coeftest(model, vcov. = cl_vcov)
  } else {cl_vcov}
  
}

#Cluster bootstrap standard errors (following Cameron, Gelbach & Miller (2008))
boot_errors = function(model, clusters, vcov = T) {
  require(multiwayvcov)
  require(lmtest)
  require(parallel)
  options(boot.ncpus = detectCores())
  cl_vcov = cluster.boot(model, clusters, 
                         boot_type = 'wild', wild_type = 'rademacher', 
                         R = 1000, parallel = T)
  
  if (!vcov) {
    coeftest(model, vcov. = cl_vcov)
  } else {cl_vcov}
}

#Write latex tables with notes:
write_latex = function(star_table, note_text, out_path) {
  table_note = note_text %>% 
    str_replace_all("\n", " ") %>%
    paste0("\\parbox[t]{\\wd\\tablebox}{", ., "}")
  out_table = append(star_table, 
                     "\\begin{lrbox}{\\tablebox}",
                     str_detect(star_table, "[\\\\]label\\{[0-9A-Za-z:_]+\\}") %>% which) %>%
    append(., 
           c("\\end{lrbox}","\\usebox{\\tablebox}\\\\[1ex]",table_note),
           str_detect(., "[\\\\]end\\{tabular\\}") %>% which)
  cat(out_table, sep = '\n', file = out_path)
}


#Simulate predicted probabilities, multinomial logit
sim_mn = function(model, newdata, sims = 1000) {
  #mnl probabilities
  mnl_p = function(logits) {
    out = vector('numeric', length = length(logits))
    for (i in seq_along(logits)) out[i] = exp(logits[i]) / sum(exp(logits))
    out
  }
  #Calculate
  mu = coef(model) %>% t %>% as.vector
  Sigma = vcov(model)
  sim_betas = mvrnorm(n = sims, mu = mu, Sigma = Sigma)
  newdata_mat = as.matrix(newdata)
  X = model.matrix(~ newdata_mat) 
  
  sim_fd = matrix(0, ncol = sims, nrow = nrow(coef(model)) + 1)
  rownames(sim_fd) = model$lab
  for (i in 1:sims) {
    beta_i = sim_betas[i,] %>% matrix(ncol = ncol(coef(model)), byrow = T)
    beta_i = rbind(rep(0, ncol(coef(model))), beta_i)
    logits = apply(X, 1, function(x) x %*% t(beta_i))
    p_hat = apply(logits, 2, mnl_p)
    fd = p_hat[,2] - p_hat[,1]
    sim_fd[,i] = fd
  }
  return(sim_fd)
}
